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summary 

A three-dimensional boimdary-layer code has been developed for particular appli- ' 
cation to realistic hjnpersonic aircraft, but it is very general and can be applied to a wide 
variety of boundary-layer flows. Laminar, transitional, and fully turbulent flows *6f com- 
pressible, reacting gases can be efficiently calculated by use of the code. A body- 
oriented orthogonal coordinate system is used for the calculation and the user has com- ' 
plete freedom in specif 3 ^ng the coordinate system within the restrictions that one. ’’ 
coordinate must be normal to the surface and the three coordinates must be mutually / : 

orthogonal. / 

The boundary -layer equations are discretized and integrated step by step. The 
Integration is fully implicit in the streamwise direction, a condition that is especially 
important for realistic configurations since it enables one to calculate flows having 
cross -flow attachment and detachment lines off the pitch plane. The code is restricted 
to those hows which are adequately represented by the boundary-layer equations; the 
analysis must be extended to adequately describe flows in which cross -flow diffusion 
"effects are important near cross-flow attachment lines. The numerical algorithm.- 
includes splined functions for dependent variables between nodes to minimize the number 
of nodes normal to the surface. This condition results in an extremely efficient solution 
procedure for reacting boundary layers. Finally, the code includes the capability to 
account. for surface normal entropy gradient effects. 

INTRODUCTION 

The design of aircraft and aerospace vehicles requires consideration of the . : r 

" Inviscid/vlscous flow field over the vehicle. Historically, the aerothermal environment ; 
of such aircraft in realistic flight conditions has been predicted by synthesizing empirical 
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data from ground -level test facilities with simplified analyses to correlate the test resets 
with the anticipated flight environment. This approach depends on the development of , 
scaling laws between ground and flight conditions. As flight regimes are further extended 
in altitude and speed to those appropriate to space shuttle and h 3 npersonic research air- 
craft, the range of applicability of scaling laws and correlations relating results from 
existing test facilities to flight conditions becomes more uncertain. The cost of con- 
structing test facilities which will yield data of sufficient quality and range to duplicate 
realistic vehicles and trajectories is prohibitive. 

Parallel with the increased flight range typified by spacecraft, space shuttle, and 
hypersonic research aircraft has been a rapid increase in the ability to obtain numerically 
exact solutions to the complete gas dynamic equations for the inviscid/viscous flow field 
about such configurations. Currently, under contract NAS 1-11525 (ref. 1), The National 
Aeronautics and.Space Administration is supporting the development of a computer code 
which will obtain the inviscid flow field solution. To predict adequately the viscous drag 
and he^t flux distributions, one must obtain the solution of the three-dimensional boundary- 
layer, equations describing the viscous layer adjacent to the vehicle surface. In this p£q>er 
recent efforts which have culminated in the development of a computer code for solving 
general three-dimensional laminar and turbulent chemically reacting botmdary layers are 
summarized. 

The numerical approach used is based on the results of two pilot codes for calcu- 
lating three-dimensional boundary layers on pointed cones at angle of attack as well as 
broad experience in developing and applying boundary -layer codes in two (dimensions. 

Key elements in the approach are: 

(a) l^lined functions to minimize nodes through the boundary-layer thickness 
(b) Finite differences used for cross -flow derivatives - - - - ; -- * 

(c) Physical variables, simplifying analyses, and changes in turbulent transport 

prc^erties model 

(d) Fully implicit solution procedure; thereby the convenient calculation of flows 

' having off -pitch-plane attachment and detachment lines is permitted 

(e) Relative ease in incorporating a wide variety of turbulence models and boundary 

conditions, including entropy layer, specified surface temperature or heat flux, 
surface catalysis, mass addition, coupling with surface ablation calculations, 
etc. 

(f) Very reasonable computing time 

The result is a code which can accurately describe laminar boundary-layer profiles with 
7 to 10 n(xles and turbulent profiles with 12 to 15 nodes through the boundary layer. The 
number of nodes through the boundary layer is important for reacting boundary layers 
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because the time required to evaluate the chemical state dominates the time spent in “ 
matrix inversion for macroscopic quantities u, v, w, and T where T is the tem- 
perature. Hence, the time implication of the number of nodes is more serious than it ’ 
would be for homogeneous boundary layers. ' Conventional implicit finite -difference 
approaches are less complicated anal57tically but require significantly more nodes in the 
plane^riormal to the wall boundary than do splined fimction procedures, simply because' 
the latter relate not only the dependent variables but also their derivatives with the nodal 
variation of the independent variables. , - < . 

Thihee -dimensional boundary layers exhibit certain physical and iriathematical 
characteristics which have important implications on the selection of a numerical solu- 
tion technique. Pairticularly important to this effort are physical and mathematical' ’ ‘ 
modeling in the vicinity of outflow and inflow lines, lines where the surface streamlines 
diverge or converge. The windward side pitch plane line of a pointed cone at incidence' ' 
is an example of an outflow line, whereas the leeward side pitch plane meridian at small 
incidence typifies an inflow line. These lines are often found off the pitch plane oii both ' 
the windward and leeward sides of configurations such as the space shuttle. They are ■ • 
important physically because they are often regions of local maximum or liiinimum heat 
flux rates. Ma.thematically, they are extremely important because of a numerical ' ' 

approach often taken in solving three-dimensional boundary-layer flows, namely, integrat- 
ing with cross flow. 

Consider the pointed cone at incidence. Because the flow direction is strictly ' 
away from the windward pitch plane (except right on the pitch plane), one expects arid 
finds that by assuming small circumferential derivatives, the solution of the boundary 
layer in the pitch plane can be determined completely independent of the boundary -layer 
solution off the pitch plane. ^ Assume for the moment that one has the entire flow field 
solution (inviscid and boundary layer) at some axial station Xq and wishes to advance 
to xi = Xq + Ax where the inviscid solution is known (Ax is the mesh point spacing). 

As noted above, the windward pitch plane solution can be obtained directly. Then, by 
recalling that the initial value problem is hyperbolic-like with respect to the circumfer- 
ential coordinate and parabolic with respect to the surface normal, the boundary -layer 
solution can be obtained by marching around the cone from the windward to tte leeward 
side. , 

In contradistinction is the class of fully implicit solution procedures in which the 
solution at Xj^ = Xq + Ax is obtained simultaneously at all points. In the fully implicit 
procedure each point at Xj is influenced by and influences every other point. The Influ- 


^Of course, one needs information about the behavior of the edge inviscid flow in the 
vicinity of the pitch plane and of the pressure, temperature, velocity components, etc. 

(that isVhot only Pg, Te, Ue, etc., but also 3pg/9<^, etc.). ' 
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ence will be lai^e or negligibly small according to the actual physical influences, pro- 
vided that the numerical scheme is accurate. 

The circumferential marching methods require that one integrate the system in the 
direction of positive cross -flow velocity in order to obey the laws of propagation of sig- 
nals. Marching against the cross flow results in instability. On a cone at very small 
incidence this is no problem, for the cross flow, is always from windward to leeward. At 
higher incidence the cross -flow inflow line moves off the leeward meridian, and it is 
necessary to march from both the windward ^d the leeward mer idians toward the- cross - 
flow inflow line. Furthermore, the marching methods can have stability limits on the 
streamwise integration step size, and the allowable step size decreases linearly with the 
meridional spacing. Fully implicit schemes usually have no such integration step size 
limit, the only limit being one of accuracy. 

Having used the pointed cone problem to clarify the situation, consider now a real- 
istic configuration such as the space shuttle. Typical outflow and inflow lines on the 
windward side are shown in sketch (a). Obviously, the solution along the pitch plane out- 



Sketch (a). - Schematic of inflow and outflow lines. 


flow line can be obtained independently of the rest of the boundary layer by assuming small 

cross -flow derivatives. This is not the cas^for the wing leading edge because there are 

no symmetry conditions there which could be used to simplify the three-dimensional 
boundary -layer equations. In fact, without solving the complete problem, one does not 
know the locus of the outflow line there.- Swept cylinder theory could be used to obtain 
an approximate solution along the wing leading edge, but this would be inconsistent with 
the spirit of an exact boundary -layer calculation. This is particularly important since 
the wing leading edge is a region of locally very high heat fluxes. The advantage of the 
fully implicit approach is ^parent, for the exact solution including the wing leading -edge 
region can be obtained without resort to approximations. 


’’ There is another very important inflow line consideration. It is well known that 
there appear to be no physically meaningful similar solutions to the full three-dimensional 
boxmdary-layer equations on the leeside of a pointed' cone at angle of attack a for values 
of a/0c certain ranges. (See refs. 2 and 3.) (0(. is the half -core angle.) One 
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expects Qomparable results to hold for inflow lines on realistic configurations.' Lin and • 
Rubin (ref. 3) suggest extending the three-dimensional boundary-layer equations to 
account for cross -flow diffusion, and they report results that show excellent agreement 
with experiment for values of a/dc which they could not obtain numerical solutions 
on the leeside by using the standard three-dimensional boundary -layer equations. At 
even higher values of oi/9q, however, they could obtain solutions with the three- 
dimensional boundary-layer equations (without cross -flow diffusion), but these solutions 
were in poor agreement with the experimental heat flux. Calculations using the extentted 
equations agreed well with the ejq>erimental data. 

ANALYSIS . 

Geometry and Coordinate Considerations 

The selection of a coordinate system to describe the development of a boundary 
layer over a three-dimensional surface is of major importance and can significantly 
affect the eventual usefulness and applicability of the final computer code. The numeri- 
cal algorithm utilized herein required the use of Taylor series spline functions (in the 
normal direction) and a Ne'wton-Raphson iteration technique to obtain the final solution 
vectors. Both of these concepts are discussed in more detail later; however, it is useful 
to emphasize that such a solution procedure requires one to generate many correction 
coefficients for each individual term in the equation set. As a result, it becomes 
extremely desirable to keep the number of terms in the equation set to a minimvun. Con- 
sequently, and because the flows to be computed will be generally highly nonsimilar, it 
was elected to retain both dependent and independent variables in their primitive forms 
(except for nodal -point stretching in the direction normal to the waH). Utilization of., , 
stream functions and similarity variables would result in significant increases in the 
number of terms in the equation set. In addition, such transformations can result in 
major difficulties as far as future code modifications and/or changes are concerned, 
especially in updating turbulence models. 

The boundary -layer approximation is impractical to implement in any other than 
surface -oriented coordinate systems. Once this is recognized, the only remaining ques- 
tion is one of orthogonal compared with nonorthogonal coordinate systems, ^gardless 
of which system is chosen, one constraint should be recognized, namely, that one coor- 
dinate direction must always be normal to the surface at all times. As a result, if one 
selects an orthogonal coordinate systemj only one other coordinate direction is indepen- 
dent, as the third will automatically be determined from the orthogonality condition. For 
example, if the streamwise coordinate is selected (for example, to be coincident to the 
body cross sections)^ the nodal -point distribution in the circumferential direction is autor 
matically determined. However, if the nodal point distribution in the circumferential 
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direction is specified, then constant streamwise coordinate stations will be, in general, 
skewed surfaces and will not lie along the body cross-section planes. On the other hand, , 
if a nonorth(^onal coordinate system is chosen, one may select both streamwise and cir- 
cumferential nodal point distributions. To decide whether to stay orthc^onal or nonor- 
thogonal requires one to investigate the form of the equations in both instances. Consi- 
der the metric tensor gy in both coordinates: 

For orthogonal systems: 



For nonorthogonal systems: 


^ hj^' hjh 2 cos 0. o\ 

gy = hjh2 cos 9 h2^ 0 

\ 0 0 1 / 

wheire hj and h2 are the scale factors. 

The nonorthogonal metric tensor results in only two additional terms in the conti- 
nuity equation beyond its orthogonal counterpart. In the momentum and energy equations* 
however, there is a threefold increase in the number of terms. By assuming that each 
terna generates, on the average, three correction coefficients, one can easily assess that 
the nonorthc^onal coordinate system results in a nominal tenfold increase in the number 
-of correction coefficients,^ - ; 

Although the desirability of choosing both surface coordinate point distributions is 
enticing, the additional complexity added to the analysis is not felt to be warranted. As a 
result, amorthogonal coordinate system was chosen; It was elected to have the flexibil- 
ity of choosing a circumferential nodal -point distribution since it is necessary to be able , 
to concentrate nodal points in regions of large cross -flow gradients (for example, wing, 
leading edges). Since the other coordinate direction is normal to the wall, the third 
• (streamwise) coordinate direction is automatically determined from the cross product. 

As a result, constant streamwise surfaces will not lie along body cross sections. This 

^This result is only true if one considers cross -flow diffusion terms within the 
equation set. For no cross -flow diffusion, the increase is minimal. Although the equa- 
tions presented herein do not include these additional terms, it is expected that their 
inclusion in the near future will be necessary in order to handle the cross-flow separa- 
tion problem and for this reason the argument has been presented. . 
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condition is inconvenient in that more handling and interpolation of e^e conditions will 
be required to obtain a consistent set of edge boundary conditions. 

It should be emphasized that the particular orthogonal system described is used in 
automatically coupling the boundary -layer code to the three-dimensional iriviscid super- 
sonic flow field code developed by Marconi et al. (ref. 1) for application to space shuttle 
problems. The boundary -layer code is written in general orthogonal surface coordinates, 
and one has the freedom to select any coordinate system within that class (for example, 
streamline coordinates), provided one can provide the inviscid edge data to the boundary- 
layer code. 


Dependent and Independent Variable Selection 

The governing equations are solved in primitive variables, but with stretching the 
independent variable normal to the surface for computational convenience. Transforma- 
tions to similarity variables are useful in anal 3 ^ically studying similarity solutions. For 
the numerical solution of nonsimilar solutions their usefulness is of questionable value, 
provided nodal points are sensibly selected and streamwise derivatives are accurately 
represented. Cebeci et al. report (ref. 4) that fewer nodes were necessary for calculat- 
ing turbulent boundary layers when they used a similarity transformation of the dependent 
variables, nodal spacing increasing geometrically from the wall. They would have had 
better results had they varied the nodal spacing logarithmically near the wall as the solu- 
tion varies. Their experience simply indicates that geometrically spaced nodal distribu- 
tions coupled with transformed dependent variables resulted in a nodal spacing which was 
closer to logarithmic in physical space than a straight geometric progression in physical 
space. To illustrate this condition, results of the 3-D code using physical variables and 
8 nodes agrees very well with the BLIMP code (ref. 5) predictions using similarity varia- 
bles and 7 nodes, as shown in figure 1. To minimize the number of points required, it is 
necessary to identify the edge of the boundary layer. This is done by scaling the lengths 
normal to the surface by the local boundary -layer thickness, an unknown quantity, which 
is obtained by introducing the auxiliary condition that the total enthalpy at the sec- 
ond (or third) nodal point from the edge is a fixed ratio of the local edge enthalpy (for- 
example, = 0.9 at X 3 /X 3 g = 0.5 where e denotes edge]. 

Governing Equations 

The three-dimensional boundary -layer equations are written in general, orthogonal^ 
body surface coordinates. One coordinate is normal to the local tangent plane, one is In 
the tangent plane in the general direction of the edge velocity, 3 and the third is in the tan- 
gent plane normal to the other two. As noted, coordinate stretching is performed normal 
to the surface. Within this context, the three-dimensional boundary -layer equations are 

^This is necessary if the problem is to be solved as an Initial value problem. 
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(subscript 1 Indicates the general streamwise coordinate; subscript 2, the cross -flow 
coordinate; and subscript 3, the surface normal coordinate): 

Continuity: , 

I- = 0 (1) 

Axial momentiun: 





where p.* is the pressure, n is the molecular viscosity, and 





£T 

ah 


and: € is the eddy viscosity. These equations are nondimensionalized, all dimensional 
quantities being divided by appropriate constant reference quantities. Important nondi- 
mensional grotqjs are listed in table I. 
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TABLE I.- IMPORTANT NONDIMENSIONAL GROUPS 


Groiq) 

symbol 


Elements 

Common name 

Re " 


P 

reference density 

Reference Reynolds number 

R 


u 

reference velocity 


Rn 

II 

£ 

reference length 

Reynolds number based on nose 



radius - 

Roo 


M 

reference viscosity . 

Free -St ream Reynolds number 

Re,“^ 

: . 



Reference free-stre^ Re 3 molds 





number 

B = 

u2 

u 

reference velocity 



h 

h 

reference enthalpy 

■' . ■■ 


Cne 

Cp 

specific heat 

Turbulent Prandtl number 

^pr,t = -r 

e 

eddy viscosity 





thermal conductivity 


f ■ 


€ 

eddy viscosity 

Turbulent Schmidt number 

Nsc,t 

pDi2 

P 

density 

- 


' ' 

Di2 

binary diffusion coefficient 

* 


Numerical Procedures 

To simplify the discussion, consider a three-dimensional Cartesian coordinate sys- 
tem x,y,z with velocity components u, v, and w. The x -coordinate corresponds to 
the direction of integration, y is locally normal to the surface, and z designates the 
cross -flow direction. The solution domain is covered by a nodal network as schematized 
in figure 2. It is assumed that the entire flow is known at some value x = X( and is to 
be determined at = X[ + Ax. For a solution at x^+j to be obtained, which satisfies 
the governing partial differential equations (PDEs) and the imposed boundary conditions, 
the functional form of the x, y, and z derivatives is specified and is substituted into 
the PDEs along with the boundary conditions. The result is a system of algebraic rela- 
tions between unknown dependent variables at the nodal points at station x = - 

Essentially all numerical solutions to the PDEs reduce to this same process. What dis- 
tinguishes one procedure from another is the functional form chosen to relate nodal values 
of dependent. variables and their derivatives, one to another,. 
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Normal to the wall boundary, the dependent variables are represented by a splined 
Taylor series between each node. Letting j denote a nodal index running from 1 at the 
wall to JMAX at the outer edge of the boundary layer results in 


2 

^+1 = Uj + u! Ay + UV ^ 




(6) 


In addition, the PDEs are integrated with reject to y between each node which insures 
that the conservation laws are satisfied exactly between nodes and greatly simplifies the 
calculation of diffusion terms by eliminating second derivatives. 


In the cross -flow direction where variations are generally much less severe than 
in the normal direction, second-order -accurate centered finite differences are used to 
represent the cross-flow derivatives; for example, 3U/8z for equal spacing is repre- 
sented as 



Uj,k+1 - Uj,k-1 
2=k+l - ^k-1 


( 7 ) 


where k is the index associated with the z coordinate. 


Axial derivatives dU/dx are represented by simple backward difference quotients 
of the form 


Uv = 




Xi+1 


( 8 ) 


where i is the index associated with the x coordinate. 

The solution is obtained fully implicitly. Thus, equations (6) and (7) are evaluated 
at station x^^i. Upon introduction of the boundary conditions ^f or example, 

^i+i JMAX k “ Ue(x,y,z), etc.^, the result is a system of nonlinear algebraic equations 
for primary (u, p, etc.) and secondary (u', p', etc.) variables at each nodal point.^ This 
system is solved by Newton -Rap hson iteration as outlined in the following paragr^h. 

Consider the system of nonlinear equations to be represented by 

Fm(U„) = 0 (9) 


where m is an equation index and Un is one of the dependent primary or secondary 
variables (Ui+ij^k, etc.j. Assume a value of » Ug and e3q>and Fm (see eq. (9)) in 
a Taylor series about Ug 

^Actually, obtaining the solution is usually enhanced if the derivative terms (u', p', 
etc.) are treated as the primary variables and the primitive variables (u, p, etc.) as the 
secondary variables. ' 
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Fm(U„) = F„(U5) + y (Uj - ug— Sa = 0 


This equation can be written in matrix form as 
A(U - U*) + F = 0 


• - ( 10 ) 


( 11 ) 


or 


where 


U = U* - A-^F 


(12) 


A-1 = 




(13) 


Thus, the solution basically reduces to tlie problem of inverting a large matrix, and 
thereby solving a system of linear algebraic equations for the unknown U. 

The matrix A is composed of a number of submatrices, illustrated in figure 3 for 
the case of five meridional (cross-flow) planes. It is block tridiagonal because cross - 
flow derivatives are represented with centered finite differences (if splined polynomials 
were used in the cross -flow direction as well, A would be dense matrix). The subma- 
trices along the diagonal are relatively dense because they include information which is 
transmitted through the boxmdary layer as well as that which is transniitted in the cross - 
flow direction. The off-diagonal submatrices are relatively sparse because they only 
contain cross -flow derivative information. 


A schematic of a typical diagonal submatrix for one meridional plane is shown in 
figure 4. Denote the submatrices by A^j and decompose the unknown vector U into 
subvectors Ui. For a given value of U denote the error in the equation by E, which is 
a vector which can be similarly decomposed into subvectors £[. Then the elements of the 
submatrix Ajj, which are denoted by amn> represent the rate of change of the mth equa- 
tion at.that node with respect to the nth dependent variable. The equations being consi- 
dered are the Taylor series equations useid in the spline fits of dependent variables, the 
governing PDEs, and certain constraint equations used in imposing boundary conditions.® 
Since the form of the Taylor series is invariant from node to node and station to station, 
the corresponding elements are constant. By taking advantage of this constancy, the sub- 
matrix can be further reduced as follows.® 


®For example, H/H^g = 0.9 at y/ye = 0-5. Additional constraints are appropri- 
ate for calculations involving vortical -layer effects. 

®Do not confuse the sub submatrix A above with the complete matrix A of fig- 
ure 3 and equation (11). 
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t 

Taylor 

series 

( 

or 


Submatrix reducticm 


A 

b' 

Ul 


'o' 

C 

D 

.^2 


E 2 


U 2 = (d - CA-1b)"^E2 


and 

Ui = A-^BU2 


(14) 


(15) 


(16) 


Return now to the complete system equations (9). Since A is block tridiagonal, 
it can be reduced by Gaussian reduction. In this reduction, it is necessary to invert a 
number of submatrices. Because of the properties of the submatrices (see discussion 
of eqs. (14) to (16)), it is only necessary to invert matrices whose size is 

4 X (Number of nodes) +3 

♦ f . (r 

Equations Constraint ' 

per node equations 

Thus, even though the complete matrix A is very large, the maximum matrix size 
which must be inverted is quite manageable, being typically about 35 x 35 for laminar 
flow and 63 x 63 for turbulent flow. 


Stagnation -Point Solution 


It is assumed that the flow is similar at the stagnation point. By placing the coor- 
dinate system origin at the stagnation point and orienting it near the stagnation point, the_ 
edge velocities can be expressed in the form"^ 



W. 


= 

£ 


(18) 


where U and W are constants supplied by the inviscid solution. In the vicinity of the 
stagnation point 


du _ ji_ 

dxj xj 

■ 

dw _ 
dxi 


(19) 


"^These expressions are for Cartesian coordinates. 
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and it is noted that 


ii.=o 

where g is the static enthalpy. In polar coordinates, 
Ug = U cos2 0 + W sin2 <p ^ 


( 20 ) 


( 21 ) 


We = (W - U)cos 0 sin 0 ^ (22) 

Since the cross -flow derivatives are well behaved at the stagnation point, the stagnation - 
point solution can be numerically generated by representing the streamwise derivatives 
as 

^ = diUje - d2U£_i (23) 


where 


dj = -d2 



(24) 


By taking 



d2 = 0 J 


(25) 


for and and di = do = 0 for the solution can be numerically generated 

dxj dxj * * dxj 

very near the stagnation point. The predicted velocity profiles for the stagnation -point 
flow on a 4:1 prolate spheroid are compared with Howarth's analjdical solution (ref. 6) in 
figure 5. 


Entropy-Layer Effects 

It is usually relatively simple to account for entropy-layer effects in two- 
dimensional or axisymmetric hypersonic flow over a bl\mt body. The usual procedure is 
to relate the mass flux in the boundary layer with the position on the bow shock wave 
encompassing the same amount of mass flux, as illustrated in sketch (b). That is, y^ 
is determined so that the mass flux between AA' equals that across the boundary layer 
BB'. The local edge conditions at B are then determined by isentropically expanding 
the streamline from behind the shock at A to the local surface pressure at B'.. It is 
relatively easy to account for the effects of normal entropy gradient as well. Naturally, 
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the solution at BB' must be obtained iteratively in order to simultaneously satisfy the 
mass balance as indicated. 



Sketch (b).- Schematic of 2-D or axisymmetric mass flux balance used 
in accounting for entropy layer effects. 


In three dimensions the problem is somewhat more complicated, but not tremen- 
dously so. ^ain, for ease of exposition, revert to simple Cartesian coordinates. The 
idea is to match the viscous and inviscid solutions in some common region of overly. 
To this end, define a stream surface S in the inviscid region, where the boundary- 
layer equations are still valid. Conservation of mass gives 
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then 


3x 


+ 


3AFz 

dZ 




( 29 ) 


Note that F - = Constant outside the boundary layer (that is, in the overlap region). 

Thus AFx and AF 2 are independent of the location of S and are functions of x 
and z only so long as S is chosen to be in the overlap region. Equation (29) and one 
of equations (28) are the relations necessary to account for entropy -layer effects in three 
dimensions. They are implemented in the . code as follows. 

From the inviscid solution the values of the AF^ . and AFz functions at each 
(x,z) location are determined; next, the edge conditions are determined as functions of 
AFx And AFz, namely 

Ue = fi(x,z,Fx^ = Fx + AFx) (30) 

(31) 

(32) ' 

AFz = Fz - Fz' = Fz - f2(x,z,Fx,AFx) (33) 


Fz^ - l2(^>A,Fx^ = Fx + AFx) 


etc. The two new equations at each meridional plane are introduced 


3AFv 3AF 


z _ 


ax 


az 


= (pv) 


wall 


where Fx and Fz are available within the existing B/L solution logic. Equa- 
tions (27) and (30) are the two constraint equations associated with entropy -layer effects 
previously discussed. 


Turbulence Modeling 

Mean field turbulence models are utilized in the present paper. The turbulent eddy 
viscosity model currently implemented in the code is the Bushnell -Beckwith model (ref. 7) 
as extended to three-dimensional flows by Harris. The eddy viscosity is given by 


where 



£ 

6 


( 35 ) 



and 


D wall damping function ; , 

y y -direction intermittency fxmction 

In the transition region, e is multiplied by a streamwise intermittency function. The 
turbulent Lewis and Prandtl numbers are taken to be constant. 

- - -- Introduction of- different -turbulent eddy viscosity models is a straightforward task, 

but it does require the derivation of certain terms needed in the’ Newton-Rap hson itera- 
tion (for example, be/ bp, ae/3u, etc.). 

Thermochemistry 

The code has ideal gas and chemical equilibrium options. In the latter case, the 
state relation is obtained by solving exactly for the species equilibrium relations in a 
special, optimized air equilibrium routine. Details are given in reference 5. It is a 
straightforward matter to replace the air equilibrium routine with a Mollier fit, if 
desired. 

PRELIMINARY RESULTS 

A number of test cases are currently in progress. In this section the results of 
predictions are compared with experimental heat-transfer data for two three-dimensional 
flows. 

The first is for a laminar boundary layer on a 10° half -angle cone at an angle of ' 
attack-of-49. The experimental data were obtained by Tracy (ref. 8) at" Moo = 7. 95^ and - 
Rg^oo = 1.25 X 10® per foot. The inviscid edge data for the boundary -layer calculation 
were obtained from Jones' tables (ref. 9). The predicted circumferential distribution of 
heat -transfer coefficient is compared with the experimental data in figure 6, where it is 
seen that the agreement is excellent. 

The second sample problem is for laminar flow over a 2.79-cm (1.1-in.) Rjj; 
sphere /150 cone at an angle of attack of 10°. The tests were conducted at free-stream 
Mach and Reynolds numbers of 10.6 and 4.1 x 10® per meter (1.2 x 10® per foot) (Cleary, 
ref. 10). Theoretical heat -transfer distributions are compared with experimental data 
in figvires 7 and 8. Axial distributions along different meridional rays are given in fig- 
ure 7 and the circumferential distribution at one fixed axial station is given in figure 8. 
The theoretical calculations were generated by using normal-shock entropy. In general, 
the agreement between theory and experiment is very good. Slightly higher predictions 
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on the windward side at x/Rjjr< 2 probably result from a father large axial step size 
and the resulting inaccuracies in numerically calculatlhg the local edge pressure 
gradient. 

CONCLUDING REMARKS 

The three-dimensional boundary -layer code reported herein has been developed 
for particular application to realistic h3^efs6nic aircraft, but it is very general and can 
be applied to a wide variety of boundary-layer flows of interest. Laminar, transitional, 
and fully, turbulent flows of compressible, reacting gases can be efficiently calculated by 
use of the code. The calculation is performed in a general body-oriented orthogonal 
coordinate system. The user has complete freedom in specifying the coordinate system 
within the restrictions that one coordinate must be normal to the surface and the three 
coordinates must be mutually orthogonal. 

The boundary -layer equations are discretized and integrated step by step. The 
integration is fully implicit in the streamwise direction. This is especially important 
for realistic configurations since it enables one to calculate flows having cross -flow 
attachment and detachment lines off the pitch plane. The code is restricted to those 
flows which are adequately represented by the boundary -layer equations; the analysis 
must be extended to describe adequately flows in which cross -flow diffusion effects are 
important near cross-flow attachment lines. The numerical algorithm includes spllned 
functions for dependent variables between nodes to minimize the number of nodes normal 
to the surface. This condition results in an extremely efficient solution procedure for 
reacting boundary layers. Finally, the code includes the cjq)ability to account for sur- 
face normal entropy gradient effects. 

Currently, the code is being coupled to the siq>ersonic, inviscid flow field code 
which was developed by Marconi et. al., Grumman Aerospace Corporation, under con- 
tract to the National Aeronautics and Space Administration. Upon completion of the 
invlscld/viscid code coiqpling, extensive comparisons for more complex configurations 
will be made. 
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Figure 1.- Illustration that physical variables require no more nodes 
than transformed variables. Laminar flow; M*, = 8; y = 1,4; 

Roo = 8.5 X 10®; 0Q = 10°; a = 0^; y, ratio of specific heat; 
shear stress at wall; 5*, displacement thickness; and q, heat 
transfer. 
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Figure 3.- Schematic of form of matrix A 
for the case of five meridional planes. 
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Figure 4.- Schematic of submatrix. 
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Figure 5. - Comparison of exact and approximate general 
3-D stagnation -point solution, v is kinematic viscosity.. 
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Figure 6.- Comparison of 3-D boundary -layer solution to experimental data. 
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^ NUMERICAL RESULTS (NORMAL SHOCK ENTROPY) 

o EXPERIMENTAL DATA 



Figure 7.0 Comparison of numerical results with experimental data. 
Axial heat-transfer distributions for 15° cone; Moo = 10.6; 

R« = 1.2 X 106 ; and o = 10°. 
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Figure 8.- Comparison of numerical results with experimental data. 
Circumferential heat-transfer distribution. Moo = 10.6; 

Roo = 1.2 X 10®; a = 10°; and x/r^ = 2.23. 
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